 % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
            % by vss
            if rov_average==1
                ScaleRov=vss;
            else
                ScaleRov=1;
            end
                       
            
            % ===================================
            % build the state space of ROV income
            % ===================================
            idio_dist=[];
            idio_dist_aut=[];
            idio_dist_sb=[];
            idio_dist_aut_all=[];idio_dist_aut_all2=[];
            
            for k1=0:vss
                for k2=0:vss
                    if length(incomevals)==2
                    if k1+k2==vss
                        % the vector of aggregate incomes from
                            % permutations of vss households
                            
                            S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1)];
                            idio_dist=[idio_dist; k1 k2];
                            
                            % the corresponding distribution of
                            % individual income values
                            
                            
                            % this computes the 'best' distribution of incomes
                            % of size vss-lagind
                            if vss>1
                                
                                II=zeros(1,2);
                                temp=[k1 k2];
                                
                                
                                for jj=1:lagind
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2]-II;
                                end
                                idio_dist_aut=[idio_dist_aut; [k1 k2]-II];
                                S_rov_aut=[S_rov_aut; ([k1 k2]-II)*incomevals];
                                II=zeros(1,2);
                                temp=[k1 k2];
                                
                                for jj=1:lagind-1
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2]-II;
                                end
                                idio_dist_sb=[idio_dist_sb; [k1 k2]-II];
                                S_rov_sb=[S_rov_sb; ([k1 k2 ]-II)*incomevals];
                                
                                
                            end 
                    end
                    else
                    for k3=0:vss
                        if k1+k2+k3==vss
                            % the vector of aggregate incomes from
                            % permutations of vss households
                            
                            S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                            idio_dist=[idio_dist; k1 k2 k3];
                            
                            % the corresponding distribution of
                            % individual income values
                            
                            
                            % this computes the 'best' distribution of incomes
                            % of size vss-lagind
                            if vss>1
                                
                                II=zeros(1,3);
                                temp=[k1 k2 k3];
                                
                                
                                for jj=1:lagind
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2 k3]-II;
                                end
                                idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                                S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                                II=zeros(1,3);
                                temp=[k1 k2 k3];
                                
                                for jj=1:lagind-1
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2 k3]-II;
                                end
                                idio_dist_sb=[idio_dist_sb; [k1 k2 k3]-II];
                                S_rov_sb=[S_rov_sb; ([k1 k2 k3]-II)*incomevals];
                                
                                
                            end
                        end    
                        end
                    end
                end
            end
           
            
            % resort S_rov in ascending order
            [S_rov,S_rovind]=sort(S_rov);
            idio_dist=idio_dist(S_rovind,:);
            
            if vss>1
                idio_dist_aut=idio_dist_aut(S_rovind,:);
                S_rov_aut=S_rov_aut(S_rovind);
                idio_dist_sb=idio_dist_sb(S_rovind,:);
                S_rov_sb=S_rov_sb(S_rovind);
            end
            
             %%%%and create S_rov_aut vectors that drops lagind individuals from
                    %%%%the rov vector. Why lagind? The biggest stable group is
                    %%%%vss+1-lagind, hence, there will be vss+1-lagind-1
                    %%%%invididuals from the rov in it, i.e. vss-lagind, hence
                    %%%%lagind individuals will not be in it.
                    if vss>1
                        for k=1:length(STABLE)
                            if STABLE(k)==1
                                for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            if length(incomevals)==2
                                                    
                                                if j1+j2==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2)
                                                    index=index+1;
                                                    temp=[j1 j2];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                                
                                            else
                                            
                                            for j3=0:vss-k
                                                
                                                if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                                                    index=index+1;
                                                    temp=[j1 j2 j3];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                                end
                                            end
                                        end
                                        end
                                    end
                                end
                            end
                            
                        end
                    
            
            % ===================================
            % simulate aggregate income transitions
            % ===================================
            tic
            T_trans=50000;
            aggincind=[];
            P_rov=[];
            
            for jjj=1:length(S_rov)
                jjj;
                inccount=[];
                agginc=[];
                for iii=find(idio_dist(jjj,:)>0)
                    
                    for kkk=1:idio_dist(jjj,iii)
                        if length(incomevals)==3
                        rng(jjj*100+iii*10+kkk,'twister');    
                        shock=rand(T_trans,1);
                        rr=(find(shock<=cumQ(iii,1)));
                        inccount(rr,kkk,iii)=1;
                        rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                        inccount(rr,kkk,iii)=2;
                        rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                        inccount(rr,kkk,iii)=3;
                        elseif length(incomevals)==2
                        rng(jjj*100+iii*10+kkk,'twister');    
                        shock=rand(T_trans,1);
                        rr=(find(shock<=cumQ(iii,1)));
                        inccount(rr,kkk,iii)=1;
                        rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                        inccount(rr,kkk,iii)=2;
                        end
                    end
                    
                    agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                end
                agginc=sum(agginc,2);
                for jjj2=1:length(S_rov)
                    qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                    aggincind(jjj,jjj2)=sum(qqq);
                end
                P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
            end
            disp('Simulation for agg income process takes')
            toc
            
            
            %%scale everything correctly
            if vss>1
                S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
                for g=vss_set(find(vss-1>=vss_set))
                            if size(STABLE,1)>0 && STABLE(g)==1
                                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                            end
                        end
            end
            
            S_rov=S_rov/ScaleRov;
            
            % ===================================
            % state space: combinations of indiv and village income
            % ===================================
            S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
            dist=[kron([1:length(incomevals)]',ones(size(idio_dist,1),1)),kron(ones(size(incomevals)),idio_dist)];
            idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
            S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
            S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
            
            if vss>1 
                        for g=vss_set(find(vss-1>=vss_set))
                            if size(STABLE,1)>0 && STABLE(g)==1 
                                for k=1:size(S_rov_aut_all,3)
                                    S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
                                end
                            end
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                    end
            
            P=kron(Q1,P_rov);
            
            Y=S(:,1)+ScaleRov*S(:,2);				% possible aggregate income outcome
            state=length(S);                    % number of states
            agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
            STATE(vss)=state;
            
            % ===================================
            % Compute set of autarky values
            % ===================================
            
            for k=1:agent
                caut(:,k)=S(:,k);
            end
            waut=zeros(state,agent);
            waut_avg_lag=zeros(state,agent);
            waut_sb=zeros(state,agent);
            % this is just the PDV of utility from consumption of
            % income caut
            
            eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
            eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
            eeewaut_rov_avg= inv(diag(ones(1,state),0)-beta.*P)*dist(:,2:end)*utility(rho,incomevals)/vss;
              % calculate the autarky value for ROV from average utility
                        % of its members in village of size vss
                        if vss>1 && Rov_sb==1
                            for j=1:state/length(incomevals)
                                % loop through individuals in the ROV
                                if lagind==1
                                    for jj=1:length(idio_dist(j,:))
                                        if idio_dist(j,jj)>0
                                            % distribution of incomes in ROV without one
                                            % individual
                                            idio_dist_lag=idio_dist(j,:);
                                            idio_dist_lag(1,jj)=idio_dist_lag(1,jj)-1;
                                            % make a temporary new 'state' of individual
                                            % and rov
                                            S_lag_temp(jj,2)=idio_dist_lag*incomevals;
                                            S_lag_temp(jj,1)=incomevals(jj);
                                            % find the corresponding state in a village of
                                            % size vss-1
                                            S_lag_2_temp=( ones(length(S_lag)/3,1)*idio_dist_lag==IDIO_DIST_ROV_KEEP(1:length(S_lag)/3,:,vss-1)) ;
                                            S_lag_2=sum(S_lag_2_temp,2)==3;
                                            S_lag_2=[S_lag_2;S_lag_2;S_lag_2];
                                            S_lag_1=(S_lag_temp(jj,1)==S_lag(:,1));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            % construct utility of individual jj
                                            if one_per_aut==1
                                                waut_temp(jj,1)=utility(rho,S_lag_temp(jj,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                            else
                                                waut_temp(jj,1)=waut_lag(S_lag_ind,1);
                                            end
                                        else
                                            waut_temp(jj,1)=0;
                                        end
                                    end
                                    
                                    % utility of individuals is average of its members
                                    waut_rov(j,1)=idio_dist(j,:)*waut_temp/sum(idio_dist(j,:));
                                    
                                    % for the second best outside options, we don't look at cases with 'holes' in the sustainable
                                    % coalitions
                                elseif lagind>1
                                    waut_rov(j,1)=100000;
                                end
                            end
                            waut_rov_sb=[waut_rov;waut_rov;waut_rov];
                        end

                        % in first iteration, or for individual deviations, waut is just consumption of income
                        if vss==1 || coaldev==0
                            for j=1:state
                                if one_per_aut==1
                                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                                    % autarky values for the rest of village are always the same - perfect
                                    % risk-sharing among the rest of village
                                    if Rov_sb==0 || vss==1
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                    else
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                    end
                                else
                                    waut(j,1)=eeewaut(j,1);
                                    % autarky values for the rest of village are always the same - perfect
                                    % risk-sharing among the rest of village
                                    if Rov_sb==0 || vss==1
                                        waut(j,2)=eeewaut(j,2);
                                    else
                                        waut(j,2)=waut_rov_sb(j,1);
                                    end
                                end
                            end
                        else
                            % this calculates autarky values for the coalitional
                            % deviations case
                            
                            
                            for j=1:state
                                % indicator for the 'autarky' state in the state
                                % space of the previous iteration S_lag
                                S_lag_2=(S_rov_aut(j)==S_lag(:,2)) ;
                                S_lag_1=(S(j,1)==S_lag(:,1));
                                S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                % for the individual, autarky value is utility of
                                % current consumption plus relevant trans
                                % probability times waut_lag, the vector of values
                                % from the next smaller insurance group
                                if one_per_aut==1
                                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                    % for rest of village, the autarky value is just
                                    % expected disc utility of current consumption
                                    if Rov_sb==0
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                    else
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                    end
                                else
                                    waut(j,1)=waut_lag(S_lag_ind,1);
                                    if Rov_sb==0
                                        waut(j,2)=eeewaut(j,2);
                                    else
                                        waut(j,2)=waut_rov_sb(j,1);
                                    end
                                end
                                
                                %%%%%Here comes the big code that cycles through all
                                %%%%%potential coalitions of people
                                for g=1:size(S_rov_aut_all,4)
                                    if size(STABLE,2)>0 && STABLE(g)==1
                                        for k=1:size(S_rov_aut_all,3)
                                            S_lag_2=(S_rov_aut_all(j,1,k,g)==S_LAG(:,2,g)) ;
                                            S_lag_1=(S(j,1)==S_LAG(:,1,g));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            
                                            agentstate=[];
                                            for kk=1:length(incomevals)
                                                if idio_dist_aut_all(j,kk,k,g)>0
                                                    agentstate=[agentstate kk];
                                                end
                                            end
                                            
                                            if size(S_lag_ind)>0
                                                for kk=1:length(agentstate) %%%note, you only need one entry per income state, not one entry per agent
                                                    
                                                    if one_per_aut==1
                                                        waut_check_ind(j,1,k,g)=utility(rho,S(j,1))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),1,g);
                                                        waut_check(j,kk,k,g)=utility(rho,incomevals(agentstate(kk)))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),2,g);
                                                    else
                                                        waut_check_ind(j,1,k,g)=WAUT_LAG(S_lag_ind,1,g);
                                                        waut_check(j,kk,k,g)=WAUT_LAG(S_lag_ind,2,g);
                                                    end
                                                    S_check(j,kk,k,g)=agentstate(kk);  
                                                    
                                                end
                                            end
                                            
                                        end
                                    end
                                end
                            end
                            
                            
                        end
                        
                        % ====================================================
                        % Introduce some penalties to achieve stability
                        % ====================================================
                        penalty=0;
                        %for village size 3
                        if Qbeta==1 & vss==2 & length(incomevals)==3 & coaldev==1
                        penalty=0.04;
                        elseif (Qbeta==2 | Qbeta==3)& vss==2 & length(incomevals)==3 & coaldev==1
                        penalty=0.03;
                        elseif (Qbeta==4 | Qbeta==5) & vss==2 & length(incomevals)==3 & coaldev==1
                        penalty=0.01;
                        elseif Qbeta==5 & vss==2 & length(incomevals)==3 & coaldev==1
                        penalty=0;
                        end
                        % for village size 4, first looking at subgroup of
                        % size 3
                        if Qbeta==3 & vss==2 & length(incomevals)==2 & coaldev==1
                        penalty=0.02;
                        elseif Qbeta==4 & vss==2 & length(incomevals)==2 & coaldev==1
                        penalty=0.0;
                        elseif Qbeta==5 & vss==2 & length(incomevals)==2 & coaldev==1
                        penalty=0.0;
                        elseif Qbeta==6 & vss==2 & length(incomevals)==2 & coaldev==1
                        penalty=0.0;
                        end
                        %for village size 4, looking at village
                        if Qbeta==3 & vss==3 & length(incomevals)==2 & coaldev==1
                        penalty=0.19;   % 0.2 does not work
                        elseif Qbeta==4 & vss==3 & length(incomevals)==2 & coaldev==1
                        penalty=0.13;  %0.1 does not work
                        elseif Qbeta==5 & vss==3 & length(incomevals)==2 & coaldev==1
                        penalty=0.01;  
                        elseif Qbeta==6 & vss==3 & length(incomevals)==2 & coaldev==1
                        penalty=0.0;  %0.1 does not work
                        end
                        
                        if Qbeta>=1 & vss>=2 & coaldev==1
                        for j=1:state
                            for g=1:2
                                if waut(j,g)>=0
                                waut(j,g)=waut(j,g)*(1-penalty);
                                waut_check(j,1,g)=waut_check(j,1,g)*(1-penalty);
                                waut_check_ind(j,1,g)=waut_check_ind(j,1,g)*(1-penalty);
                                eeewaut(j,1)=eeewaut(j,1)*(1-penalty);
                                elseif waut(j,g)<0
                                waut(j,g)=waut(j,g)/(1-penalty);
                                eeewaut(j,1)=eeewaut(j,1)/(1-penalty);
                                waut_check(j,1,g)=waut_check(j,1,g)/(1-penalty);
                                waut_check_ind(j,1,g)=waut_check_ind(j,1,g)/(1-penalty);
                                end
                            end
                        end
                        end
                        % ====================================================
                        % Compute constrained insurance
                        % ====================================================
                        gammagrid=[lambda vss/ScaleRov*(1-lambda)];
                        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
                        gammagridold=gammagrid;
                        clear gammagrid
                        gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
                        
                        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
                        n=length(gammagrid);
                        w=zeros(n,state,agent);
                        gammar=zeros(n,state,agent);
                        c=zeros(n,state,agent);
                        phi=zeros(n,state,agent);
                        
                        % 'first-best' consumption given aggregate resources and planner weights
                        for j=1:state
                            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
                        end
                        cfirstbest=c;
                        %  relative planner weights
                        for j=1:state
                            for k=1:agent
                                gammar(:,j,k)=gammagrid(:,agent+k);
                            end
                        end
                        % Compute the first best value functions without binding
                        % constraints
                        w=zeros(n,state,agent); dww=1;
                        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
                        while dww>0.00001
                            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                            dww=max(max(max(abs(w-wold))));
                            wold=w;
                        end
                        wfirstbest=w;
                        % continuation utility as function of planner weight and state of the
                        % economy
                        ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
                        ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
                        
                        % ====================================================
                        % Start the loop with the enforcement constraints
                        % ====================================================n
                        ewsb=ewfirstbest;
                        wsb=wfirstbest;
                        options=optimset('Display','off');
                        eps=0.0000001;t=0;itgap=0; GAP=2; gap=2; gaphist=gap; gapincrease=0;
                        csolold=cfirstbest;
                        if beta<0.85
                            gapcrit=gapcritset;
                            maxitgap=maxitgapset;
                        else
                            gapcrit=gapcritset*1
                            maxitgap=maxitgapset;
                        end
                        
                        %%idiosyncratic distribution in rov for all states 
                        while GAP >=gapcrit && itgap<=maxitgap
                            rr=[];
                            rr1=[];
                            rr2=[];
                            rrnone=[];
                            rrboth=[];
                            nonconv=zeros(state,3)*nan;
                            itgap=itgap+1;
                            coalnonsust=nan(state,2);
                            count=0; countgrid=[ones(n,state)];
                            for j=1:state
                                
                                
                                
                                index=[1 j];
                                % 1. constraint binding for agent 1
                                % (individual)
                                rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if vss>1 & coaldev==1
                                    rr1cd=[];
                                    for g=size(waut_check,4):-1:0 %%%check from the largest group down, 0 is an individual deviation
                                        
                                        if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                                            %%%check the largest possible coalition first:
                                            if g>0
                                                for k=1:size(waut_check,3)  
                                                    if waut_check(j,1,k,g)~=0
                                                        coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-waut_check_ind(j,1,k,g)];
                                                        for kk=1:size(waut_check,2)
                                                            if waut_check(j,kk,k,g)~=0
                                                                coalition=[coalition utility(rho,cfirstbest(rr1,j,2))+beta*ewsb(rr1,j,2)-waut_check(j,kk,k,g)];
                                                            end
                                                        end
                                                        %%%Step 1: find the grid points for which you are now trying to find a solution
                                                        
                                                        %%%%do the others want to come along?
                                                        rrcheck=coalition<0;
                                                        rrcheck=prod(rrcheck,2);
                                                        rrcheck=find(rrcheck>0);
                                                        if size(rrcheck>0)
                                                            %%%%you need to do the solutions here
                                                            bind=1;
                                                           
                                                            coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            rrchecksolve=coalitionsolve<0;
                                                            rrchecksolve=prod(rrchecksolve,2);
                                                            rrchecksolve=find(rrchecksolve>0);
                                                            
                                                            rr=rrchecksolve(end)+1;
                                                            if rr>length(gammagrid)
                                                                coalnonsust(j,1)=-2;
                                                                output=[nan,nan,-50];
                                                            else
                                                                
                                                                clear x0
                                                                
                                                                x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                                %T vectors of consumption and utilty that meet
                                                                %participation constraints
                                                                [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[waut_check_ind(:,1,k,g) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                                                
                                                                %T 29_Mar This keeps track of nonconvergence, but
                                                                %allows the programme to continue - sometimes it
                                                                %converges in a subsequent iteration
                                                            end
                                                            if output(3) <=0
                                                                nonconv(j,1)=output(3);
                                                                csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                                csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            else
                                                                csol(rr1(rrcheck),j,1)=x(1);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                                csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            end
                                                        end
                                                    end
                                                end
                                                
                                                
                                                rr1(rrcheck)=[];  %%%these grid points have been dealt with
                                                rrnoneadd=rr1;
                                            elseif g==0 %%%only an individual deviation
                                                coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-eeewaut(j,1)];
                                                
                                                %%%%do the others want to come along?
                                                rrcheck=coalition<0;
                                                rrcheck=prod(rrcheck,2);
                                                rrcheck=find(rrcheck>0);
                                                if size(rrcheck>0)
                                                    
                                                    bind=1;
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];
                                                    %%%%do the others want to come along?
                                                    rrchecksolve=coalitionsolve<0;
                                                    rrchecksolve=prod(rrchecksolve,2);
                                                    rrchecksolve=find(rrchecksolve>0);
                                                    
                                                    
                                                    rr=rrchecksolve(end)+1;
                                                    
                                                    if rr>length(gammagrid)
                                                        coalnonsust(j,1)=-2;
                                                        output=[nan,nan,-50];
                                                    else
                                                        clear x0
                                                        x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                        %T vectors of consumption and utilty that meet
                                                        %participation constraints
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                                        
                                                        %T 29_Mar This keeps track of nonconvergence, but
                                                        %allows the programme to continue - sometimes it
                                                        %converges in a subsequent iteration
                                                    end
                                                    if output(3) <=0
                                                        nonconv(j,1)=output(3);
                                                        csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                        csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    else
                                                        csol(rr1(rrcheck),j,1)=x(1);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                        csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    end
                                                end
                                                rr1(rrcheck)=[];  %%%because these grid points have been dealt with
                                                rrnoneadd=rr1;
                                            end
                                        end
                                    end
                                else
                                    %%%%in first group size iteration, i.e. for vss=1
                                    rr1cd=rr1;
                                    rr1id=[];
                                    rrnoneadd=[];
                                end

                                
                                if size(rr1cd)>0
                                    bind=1;
                                   
                                    rr=rr1cd(end)+1;
                                   
                                    if isempty(rr)==1
                                        coalnonsust(j,1)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        if size(rr1cd)>0
                                            clear x0
                                            x0(1) = cfirstbest(rr1cd(1),j,1);  %%redundant because we no longer need an initial guess
                                            %T vectors of consumption and utilty that meet
                                            %participation constraints
                                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                            
                                            %T 29_Mar This keeps track of nonconvergence, but
                                            %allows the programme to continue - sometimes it
                                            %converges in a subsequent iteration
                                        end
                                        
                                    end
                                    if output(3) <=0
                                        nonconv(j,1)=output(3);
                                        csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                        csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                        gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                        wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    else
                                        csol(rr1cd,j,1)=x(1);
                                        gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=output(1);
                                        wsbnewsol(rr1cd,j,2)=output(2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    end
                                end
                                %%%%this becomes relevant for vss>1, these are all the left-overs that are not dealt with by joint deviations.
                                if size(rrnoneadd)>0
                                    
                                    csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                                    gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                                    phisol(rrnoneadd,j,1:agent)=0;
                                    
                                    wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                                    wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                                end
                                
                                % 2. Constraint binding for agent 2 (RoV),
                                rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                                
                                if size(rr2)>0
                                    bind=2; %%i.e. agent 2 has the binding constraint
                                    rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                                    if isempty(rr)==1
                                        coalnonsust(j,2)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        clear x0
                                        x0(1) = cfirstbest(rr(end),j,1); %%redundant because we no longer need an initial guess
                                       
                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                        %%%%check that agent 1 does not have a binding
                                        %%%%constraint at this new solution, and note that
                                        %%%%the rov must never have less than waut(j,2),
                                        %%%%so if individual has a binding constraint that
                                        %%%%means rov would have to decrease its
                                        %%%%consumption, not possible, so you do not need
                                        %%%%to solve, just check!
                                    end
                                    % output(3) equals 1 if constraints
                                    % are satisfied, but -2 if not
                                    if output(3)<=0
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        %%%%first check that agent 1 (in conjunction with the rest of the village)
                                        %%%%truly does not have binding constraints
                                        %%%%at the proposed solution
                                        if vss>1  & coaldev==1 % only need to check this once coalitional deviations become relevant
                                            for g=size(waut_check,4):-1:1 %%%check from the largest group down, 0 is an individual deviation
                                                
                                                if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                                                    %%%check the largest possible coalition first:
                                                    for k=1:size(waut_check,3)  
                                                        if waut_check(j,1,k,g)~=0
                                                            coalition=[output(1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalition=[coalition output(2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            %%%Step 1: find the grid points for which you are now trying to find a solution
                                                            %%%%do the others want to come along?
                                                            rrcheck=coalition<0;
                                                            rrcheck=prod(rrcheck,2);
                                                            rrcheck=find(rrcheck>0);
                                                            if size(rrcheck)>0
                                                                output(3)=-2;
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                        
                                        if output(3)>0
                                            csol(rr2,j,1)=x(1);
                                            gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                            csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                            wsbnewsol(rr2,j,1)=output(1);
                                            wsbnewsol(rr2,j,2)=output(2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        else
                                            nonconv(j,2)=output(3)*100;
                                            csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                            gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                            csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                            wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                            wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        end
                                    end
                                end
                                
                                % 3. Constraint all slack
                                rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if size(rrnone)>0
                                    csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                                    gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                                    phisol(rrnone,j,1:agent)=0;
                                    
                                    wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                                    wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                                end
                            end
                            gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                            
                            gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                            gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                            gap=max(gap1,gap2);
                            gaphist(itgap+1)=gap;
                            GAP=max(gaphist(max(1,length(gaphist)-convlength):end));
                            
                            disp('gap,gapold')
                            disp([gap,gapold])
                            nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                            Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                            
                           
                            if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                                itgap
                                disp('Autarky or coalition not sustainable')
                                [nonconv,S]
                                for j=1:state
                                    csol(:,j,1)=caut(j,1);
                                    csol(:,j,2)=caut(j,2);
                                end
                                break
                            end
                            csolold=csol;
                            wsb=wsbnewsol;
                            ewsb(:,:,1)=(P*wsb(:,:,1)')';
                            ewsb(:,:,2)=(P*wsb(:,:,2)')';
                        end
                        CRITGAP(vss)=gaphist(itgap-1);
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        if gapincrease==1
                            nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
                        end
                        
                        % =============================================================
                        % calculate values of outside option for coalitional deviations
                        % =============================================================
                        
                        if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                            % if vss is unsustainable, the outside option for vss+1 remains that
                            % for vss
                            waut_lag(:,1)=waut_lag(:,1);
                            waut_lag(:,2)=waut_lag(:,2);
                            P_lag=P_lag;
                            S_lag=S_lag;
                            % if the current insurance group is not sustainable, it
                            % cannot be the outside option for the next bigger group,
                            % so increase the 'lagind' indicator by one
                            lagind=lagind+1;
                            gap=-1;
                            STABLE(vss)=0;
                            STABLECHECK=[1 STABLE];
                            % ... or if there is no previous, use just conditional exp of
                            % autarky consumption as calculated above
                        elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                            waut_lag(:,1)=eeewaut(:,1);
                            waut_lag(:,2)=eeewaut(:,2);
                            P_lag=P;
                            S_lag=S;
                            lagind=1;
                            gap=-1;
                            
                            
                            S_rov_KEEP(1:state/length(incomevals),vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE];
                            IDIO_DIST_ROV_KEEP(1:state/length(incomevals),1:length(incomevals),vss)=idio_dist;
                        else 
                            lagind=1;
                            % when ins group of size vss is sustainable, its value
                            % becomes outside option for next bigger group size
                            for j=1:state
                                % outside option for next larger insurance group is
                                % value with equal weights vss ...
                                eewaut(j,1)=wsb((n+1)/2,j,1);
                                eewaut(j,2)=wsb((n+1)/2,j,2);
                            end
                            waut_lag=eewaut;
                            S_lag=S;
                            P_lag=P;
                            S_rov_KEEP(1:state/length(incomevals),vss)=S_rov;
                                                        
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            WAUT_LAG(1:state,1:2,vss)=waut_lag;
                            P_LAG(1:state,1:state,vss)=P_lag;
                            S_LAG(1:state,:,vss)=S_lag;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE]; %%%this includes stability for the individual as well, makes the loop above easier
                            
                            IDIO_DIST_ROV_KEEP(1:state/length(incomevals),1:length(incomevals),vss)=idio_dist;
                            
                        end
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
